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ABSTRACT 

This papers explores the self similar solutions of the Vlasov-Poisson system and their 
relation to the gravitational collapse of dynamically cold systems. Analytic solutions 
are derived for power law potentials in one dimension, and extensions of these so- 
lutions in three dimensions are proposed. Next the self similarity of the collapse of 
cold dynamical systems is investigated numerically. The fold system in phase space 
is consistent with analytic self similar solutions, the solutions present all the proper 
self-similar scaling. An additional point is the appearance of an law at the center 
of the system for initial conditions with power law index larger than — ~ (the Bin- 

ney conjecture). It is found that the first appearance of the law corresponds to 
the formation of a singularity very close to the center. Finally the general properties 
of self similar multi dimensional solutions near equilibrium are investigated. Smooth 
and continuous self similar solutions have power law behavior at equilibrium. However 
cold initial conditions result in discontinuous phase space solutions, and the smoothed 
phase space density looses its auto similar properties. This problem is easily solved by 
observing that the probability distribution of the phase space density P is identical 
except for scaling parameters to the probability distribution of the smoothed phase 
space density Ps- As a consequence Ps inherits the self similar properties of P. This 
particular property is at the origin of the universal power law observed in numeri- 
cal simulation for The self similar properties of P$ implies that other quantities 
should have also an universal power law behavior with predictable exponents. This 
hypothesis is tested using a numerical model of the phase space density of cold dark 
matter haloes, an excellent agreement is obtained. 
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1 INTRODUCTION. 

The formation of dark matter haloes in cosmology has been 
investigated in great details and impressive accuracy by nu- 
merical simulations. An important discovery was the exis- 
tence of the universal NFW density profile (Navarro, Frenk, 
& White, 1997). Later Taylor & Navarro (2001) discovered 
another intriguing property, the so-called pseudo-phase den- 
sity, -Pg where p is the projected density and a is the velocity 
dispersion, is a universal power law as a function of radius. 
The exponent of this power law corresponds exactly to the 
exponent obtained for the self-similar solutions developped 
by Fillmore & Goldreich (1984) in the special case of the 
secondary infall model studied by Bertschinger (1985). The 
initial purely radial Bertschinger model was recently refined 
by Zukin & Bertschinger (2010a, 2010b) by introducing tidal 
torque. In this new approach Zukin & Bertschinger con- 
firmed the universality of the pseudo phase space density at 
large and intermediate scale but questioned the universality 
at small scale. Very high resolution numerical simulations 
also raise questions about the universality of the pseudo 



phase space density at small scale (Stadel etal 2009). An 
interesting addition to this debate is the work of Nusser A. 
(2001) who show that the introduction of momenta does not 
modify the mass profile. Recent investigations confirmed 
the universal power law behavior of the pseudo phase space 
density, and improved the numerical accuracy by taking into 
account the ellipticity of the haloes and the substructures. 
As an example Ludlow etal (2010) demonstrated that the 
residual to the fit of a power law to the pseudo phase space 
density is typically about 10 to 20% and does not exceed 
30%. The power law behavior of the pseudo phase space 
density is now observed with impressive accuracy over 3 
decades in radius. Note that the universality of power laws 
for the pseudo-phase density is stronger than the univer- 
sality of the NFW profile (see Vogelsberger, Mohayaee, & 
White 2011). For a general review of the universal and not so 
universal properties of dark matter haloes, see Navarro etal 
(2010). The relation with the Bertschinger solution and the 
strong power law behavior of the pseudo phase space den- 
sity seems to indicate a relation to the self similar solutions 
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of the Vlasov Poisson system. This paper will explore the 
details of the relation between the collapse of a dynamically 
cold system in a cosmological context and self similarity. Be- 
fore tackling the general problem of the multi dimensional 
solutions, some attention will be given to a simpler problem, 
the one dimensional cold collapse. In particular it is inter- 
esting to consider and study another universal behavior, the 
convergence of the one dimensional cold collapse to a uni- 
versal power law density with exponent — | for a large class 
of initial conditions (Binney 2004). 



2 THE VLASOV-POISSON SYSTEM. 

The Vlasov-Poisson system for a space space density 
/(x, v,t), with general space coordinates x velocity v, time 
t, potential <?!>(x) and density p(x), is: 



21 i 9/ v _ 2ia± 

dt dx dv dx 

A(j) = q n Gp 

P = ffd"v 
with q„ a constant. 



= 



(1) 
(2) 
(3) 



3 SCALE FREE SOLUTIONS. 

The scale free solutions should be consistent with the fol- 
lowing condition: 



/(Aix, A 2 v, A 3 i) = A 4 /(x, v, t) 



(4) 



Considering the special case, Ai = 1 + edAi, and expanding 
to the first order in e, the following equation is obtained: 

|^dA 1 x+^-dA 2 v + ^dA 3 t-/dA 4 = (5) 
ax Ov Ot 

The solution to Eq [5] reads: 

/(x,v,t) = /ot- F (_^,_y 

with: 

dA 4 dAi dA 2 

ao = -7T— , Qi = -rr~, Q2 = 



(6) 



dA 3 ' 



dA 3 ' 



dA 3 



4 SCALE FREE SOLUTION OF THE VLASOV 
EQUATION. 

The self similar solutions of the Vlasov equation are found 
by inserting Eq [5] in Eq [T] We define the following scaled 
variables: X2 = t^t, V2 = j^j. Using the former notations, 
a solution of the type given in Eq[6]is a solution of the Vlasov 
equation if: 

OF OF 
(2 + na 2 )F + (1 + a 2 )- — x 2 + tt2- — v 2 
CTX 2 9v 2 



_ / OF &4> dF 

\<9x 2 <9x 2 0v 2 



(7) 



and: 



ao = — 2 — na 2 , ai = 1 + a 2 (8) 
The definition of the scaled potential reads: 



Note that an equation similar to Eq[8]was already derived 
by Lancellotti and Kiessling (2001). 



5 ONE DIMENSIONAL SOLUTIONS. 

The one dimensional solution (n = 1 in Eq(7]) corresponds 
to a two dimensional phase space described by coordinates 
(x,v), and and associated self similar coordinates (:r 2 ,u 2 ). 
The solutions of Eq [7] will be investigated for a power law 
potential, other quantities like the force and the density are 
also assumed to be power laws. 



[X2 



k\x2 



/3+2 



(10) 



The variable X2 is re-scaled so that k = i. .The following 
change of variable u = \x 2 \ n is introduced with rj = ^ + 1 
in Eq[7] It is assumed that — 2</3<0, as a consequence 
the exponents of 4> an d u are positive. But note that this 
requirement is not sufficient since the divergence of the total 
mass must be avoided. For /3 < — 1 a finite mass would imply 
a cutoff of the power law density at the center, and would 
generate a force which would not be a power law anymore. 
This additional requirement implies that — 1 < j3 < for 
this type of one dimensional solutions. The corresponding 
equation for the variable u is: 

, OG dG 
2 + Q 2 + (1 + a 2 m— — m + a 2 - — v 2 

Qu 0V2 



,'dG „, dG\ J_ 
-sgn(x 2 )?) — — v 2 - 2ku- — 

OU OV2 J 



(11) 



where sgn(a;) is the function: 



sgn(a:) = 



-1 x < 
1 x>0 



and log (F(x2, v 2 )) = G(u, v 2 ) We introduce again new vari- 
ables, (R,ip), with: 



R = \Ju 2 + v't 
cos if) - 



B ^2>0 
— | X2 < 



The new equation in these variables is: 

OH OH 
—r/2 cost/; sin'0(l + Q2 )"^" + R ( c °s 2 ^(1 + c 2 )rj2 + a 2 ) 



+77 ^ (R\ costp\)7T2 +2 + 02 = 

where 1J2 = Tj — , and H(R, ip) — G(u, v) 

An interesting case is when the potential 4>{x) is sta- 
tionary, which corresponds to, f3ai = —2 and to r\2 — 
(see Eqs [8] and [9]) . In this case Eq[l2]is reduced to: 



OR 
(12) 



„ OH q 2 OH . 
2 + a 2 + Ra 2m + T ^- w (R\cosm\) 

Eq[l3]has a general solution: 

H(R^) = -^±l [ogR+Q ( + l+^l 
012 V 012 



(13) 



cos(^)| a2 dip 



X2 = r 



(9) 



(14) 

Eq[T2]has an asymptotic solution for small R (and all values 
of rj2), the term in Rf+ 2 that multiplies the ip derivative of 
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H dominates other derivatives of H. The solution is: 



With: 



cos^l T^dtp + H^R) 



For a stationary potential, which corresponds to j3ai 
Eq[l5]is reduced to: 



(15) 

= -2, 



H{R,1>) 



(a 2 + 1)(«2 + 2) 



a 2 



i?°2 J | cos V| Q 2 dip + H 2 (R) 
(16) 

Note that developing the general solution in Eq[l4]for r\ 2 = 
0, small R, and adopting Q(s) — — (2 + a 2 ) hi(s) + Q2(s) we 
obtain: 



H(R,tP)~Q 2 (R °2)+ 

1 

■ J* I cos ip\ " 2 



OL2 



~(a 2 + 2)P°2 +Q 2 (R °2 



Provided that Q 2 (R ° 2 ) is of lower order than P a 2 the 
former expression reduces to Eq |16l Note that if Q 2 behaves 
at origin like a power law with negative exponent this 
condition will be satisfied, and that as a consequence the 
behavior of the Q 2 solution at origin will not be singular, 
which is a condition for the validity of the approximation 
at small R. As a consequence the general solution and the 
approximation at small R are consistent for 772 = 0. 



6 GEOMETRY OF THE SOLUTION. 

To study the geometry of the solution the iso-contours of 
the solution will be investigated. Considering a contour with 
phase space density Fo, and using Eq[21 the equation of the 
corresponding contour is: 



Q(r £ = 
With the following definitions: 
Q = InQ and I(tP) = (1 + a 2 ) 



R °2 Fo 



(17) 



5(V)| "3 dV, (18) 



it is easy to realize that Eq [17] corresponds to the equa- 
tion of a spiral in phase space. Since I(tp) — ktp where k 
is a constant, the equation reduces to R — s(ip), a specific 
example is a logarithmic spiral which would occur if s is 
a power law plus a constant. An important parameter of 
this spiral structure is the inter-fold distance. This parame- 
ter is straightforward to measure in numerical simulations, 
and can be easily compared to the analytic expectation. The 
variation of tp between two consecutive folds is 2n, thus if 
the position on one fold is tpo the corresponding position on 
the consecutive fold is tpi = ipa+2ir. We will assume that the 
folds corresponds to a large number of turns, consequently 

i> > 2tt, and J(i/>i) = I(ipo) + f^ | cos(#)|^# = Io + eh, 
with e < 1. Similarly we write the inter fold variation of 
Ri — Ro + edR. Inserting these expressions in Eq 1171 and 
using the original equation we obtain: 



t»2 + 2 

dR PjR a a2 F )ha 2 

R ~ 2|±2 — — 

P(R a2 F )R ° 2 +a 2 + 2 



(19) 



P{u) = 



)' (<3 _1 («)) 



(20) 



Requiring that all functions are power law's and starting 
with Q, then Eq 1201 implies that P is also a power law. If 
^£ is also a power law Eq 1191 implies P{u) oc u a 2 , and as a 
consequence we have: 



dR 
Ro 



p»2 + l 
t 



ha 2 



-Rq 2 oc Rq ' 



(21) 

F ° 2+2 + Q 2 + 2 

An interesting consequence of Eq [5T] is the prediction of 
the projected density caustics. Introducing R — x? + and 



a 2 = — 4 — 1, the position p n of n fold is: 



(22) 



The self similar solution of Bertschinger (1985) corresponds 
to P = — I which implies p n — n~o, Sikivie & Kinney (2001) 
proposed the following numerical approximation, p n ~ 
which is quite close to the prediction of Eq 1221 For more 
details about the formation of caustics and the fine grain 
structure of phase space in the cold dark matter model see, 
Mohayaee etal (2006), Afshordi, Mohayaee & Bertschinger 
(2008), Duffy & Sikivie (2008), Vogelsberger etal (2008), 
White & Vogelsberger (2009), Vogelsberger & White (2011). 
It is important to note that the thickness of a given contour 
can be derived using the same method as the inter-folds 
distance. The variation of position for a contour in Eq 1171 
corresponds to a variation of I(il>i) = I{ipo + dtp), where 
dtp corresponds to a differential in the initial conditions for 
two opposite sides of the contour. As a consequence the cal- 
culations are identical to the calculations performed for the 
inter-fold distance, and the contour thickness is proportional 
to the inter-fold distance. 



7 MULTI DIMENSIONAL SOLUTIONS. 

It is useful to define the self similar coordinates in the 6D 
space, (x 2 , v 2 ), and the corresponding modulus of distances 
r and velocities v. The solution in several dimensions is con- 
structed using the results obtained in one dimension. Note 
that in three dimension /3 > — 1 is not required anymore, to 
avoid the divergence of the mass at center /3 > —3 is suf- 
ficient. The solution is the sum of a general solution plus 
a special solution. To facilitate the calculations we intro- 
duce G(x2,v 2 ) = logF(x2,V2), and a power law potential 
<£(r) = ir /3+2 , -2 < ft < 0, EqElthen reads: 



(2 + na 2 ) + Pi+P 2 = 
Pi = (l + a 2 )H-.X2 + a2 



8G 



V 2 



(23) 



ft = -||.v a + 08 + 2)r^.x a 

considering a solution of the type: 

G(x 2 ,v 2 ) =Ko{R)+Ki(r,v r ) (24) 

Here R = yjv 2 +r?+ 2 , and Ki(r, v r ) is a function of the 
distance modulus and radial velocity v r . The special solution 
K\ is constructed using the method we developed for the 
one dimensional solution. Provided we substitute the (|x|, v) 
space with the (r, v r ) space the equations are the same, and 
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finding the solution K\ is thus a problem we have already 
solved. Assuming that K\ is a full solution or an asymptotic 
solution for small P, we introduce the solution Eq [2H in 
Eq 1231 Thus considering G — Ko + K\ the only remaining 
significant terms may come from Ko only. Since basically Ko 
is a function of energy, we have P2 = 0. Using the P variable, 
we find that Pi behaves like R^^-. Since the terms relative 
to K\ are of order of unity, the Pi term of Kq must be 
compared to unity. Provided that the function Ko behaves 
like a power law with positive exponent, the relevant Pi term 
will become negligible at small R and the solution in Eq 1241 
is an asymptotic solution at small R. 



8 PSEUDO PHASE SPACE DENSITY FOR 
THE THREE DIMENSIONAL SOLUTION. 

The power-law behavior of the pseudo phase space density 
observed in 3D numerical simulations is a strong motivation 
to study the pseudo phase space density of the solutions. 
For a comparison to the numerical simulations the quantity 
of interest is the smooth pseudo phase space density (see 
Sec. 9 for a definition of the smooth quantities). It is simple 
to construct an analytic self similar solution at equilibrium 
by using a power law model. Provided all histograms of the 
values of / are power laws, the smoothed density of / will 
also correspond to a power law. For an equilibrium model 
it is straightforward to derive the parameters of this power 
law model. However due to the geometric nature of the so- 
lution (see Sec. 6) the cold solution occupies a finite range. 
The cold solution is a series of folds, with the first fold as 
the boundary of the system. Sec. 6 shows that this outer 
boundary may be approximated with a constant radius in 
R, R = Po- As a consequence the density of the cold solu- 
tion and of its smoothed expectation must be cut outside 
some radius Po. The smooth pseudo phase density for the 
cold density is thus obtained by estimating the expectation 
for the smooth power law model with initial density r~^° 
within Po. The variable r is the modulus of x, at this point 
it is also useful to introduce the variable v, the modulus of 
v. In three dimension the pseudo phase density is thus 
the following integrals have to be estimated numerically: 

C v ° i: 2 

P = J (y +0Vdo 

j2 (;;°(^+0)v^) 

p 

With: 

v = ^2{Rl - 4>) 

= 1 A) + 6 
V ~ 2 A) + 2 

And: 

<t> = ^r Po+2 

note that in the above integrals the re-scaling of Po implies 
the re-scaling of r and the re-scaling of the integrals, thus 
the general shape of may be obtained for any value of 
Po. The results of the calculations are presented in Fig [T] 




-5 -4 -3 -2-10 1 
log(r) 



Figure 1. The general shape of the smoothed pseudo phase space 
density for cold 3D self-similar solutions. The pseudo phase space 
density was obtained by integrating the power law expectation 
within a boundary, see Sec. 8 for more explanation. The variations 
of -£3 are represented from top to bottom for 0o = — g , Po = — § > 

A> = -f. 

The pseudo phase density of the solution is consistent with 
the theoretical power law expectation near the center and 
deviates near the edge of the system. This is obviously ex- 
pected since the cut has more effect at the edge than at the 
center. An interesting consequence of the behavior of these 
solutions, is that even if the pseudo phase density deviates 
significantly from a power law near the edge of the system 
it is still an auto-similar solution in this area. Although the 
solution computed here is isotropic which is not the case 
of the dark matter haloes investigated in numerical simula- 
tions, this result suggests that even if the pseudo phase space 
density is not a power law in some area, it is still possible 
that the phase space density is self-similar in this area. 



9 ONE DIMENSIONAL COLD COLLAPSE. 

9.1 Power law regime in the central region. 

The evolution of a system of particles with small initial ve- 
locity dispersion in the two dimensional phase space under 
the action of one dimensional planar gravity will be inves- 
tigated in this section. Interestingly Binney, (2004) found 
that for such systems the density near the center is very 
close to x~2 . This power law density suggests the possibil- 
ity of a universal self-similar regime in the central region. 
Binney (2004) proposed that the x~ 5 regime is valid for a 
very large class of initial conditions, but an interesting prob- 
lem is to find for what type of initial conditions the x~ 2 is 
not obtained, and where is the limit. This problem may be 
solved by using the following simple approach, we take a 
set of cold initial conditions with power law densities and 
compute the density of the system at late times. The initial 
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conditions are sampled by using 4 x 10 particles, the phase 
space coordinates (x,v) of these particles are initiated using 
random number generators. The x positions are initiated us- 
ing a power law random generator, while a Gaussian random 
generator is used for the v coordinates. The x dimension ran- 
dom number generator is constructed by applying a power 
law transform to the output of a uniform random generator. 
The initial conditions are evolved using a leapfrog integrator 
(Birdsall & Langdon, 1985). The leapfrog is a symplectic in- 
tegrator which is particularly important for reconstructing 
the structure of phase space. The leapfrog method requires 
the evaluation of the force at the position of each particle 
which is computed by subtracting the total mass on each 
side of the particle. No smoothing is applied to the force. 
The final step is to evaluate the proper time step for the 
leapfrog integration. The time step must be small enough to 
allow a proper reconstruction of the shortest orbits in the 
system. The shortest orbits are the closest to the center of 
the system. Assuming that the typical size of the system is 
unity if we require a dynamical range in x of Do , the associ- 
ated shortest distance is xo = . In the pre-collapse regime 
the fall time of a particle To to the center depends on the ini- 
tial mass Mo inside the particle shell. In a system with initial 
power law density p = pox 130 , T = ^ = \J %j^ x ~^- 
The proper time step for leapfrog integration should be a 
fraction of To; numerical experiments show that a fraction 
of 5% is optimal. Another point is the estimation of the ini- 
tial velocity dispersion ao- Basically the resolution in the 
velocity space ao has to match the resolution in x space 
xo- The velocity at the center of a particle with initial dis- 
tance xo should be of the order of ao, which corresponds 

to a ~ V2GM x = 2y!-Q&-Xo r+1 - Note that the above 
calculations of To and ao are not fully analytic for ft < — 1, 
in this case the mass Mo was computed directly from the 
initial conditions. The results of the ID numerical simula- 
tions are presented in Fig [2] The resolution adopted corre- 
sponds to Do — 1000. For a better numerical efficiency, in 
the range — 1 < /?o < the power law exponent of the force 
/3f is computed, the force is much smoother than the den- 
sity (affected by caustics), the equivalent projected density 
exponent, is (3f — 1- For — 2 < f3o < — 1 the power law ex- 
ponent was computed directly by fitting a power law to the 
density. Interestingly the initial exponent /3o is conserved 
at later time for /3o < — s, while for larger values of /3o the 
final exponent has a constant value of — \ (see Fig [2]). As 

a consequence the x~ 2 density at the center is not univer- 
sal but does occur precisely when fto = Note that this 
effect was also noticed by Schulz etal (2012) with results 
quite similar to Fig [2] This particular behavior suggests a 
cut-off, no final exponent is larger than — | . To understand 
the cause of this cut-off it is interesting to investigate nu- 
merical simulations with initial conditions > /3o > — i and 

look for the first appearance of the x~? power law. This 
analysis is illustrated for fio = — |, the first appearance of 

the x~ 2 projected density occurs when the the first fold 
starts to form or equivalently just after the first crossing of 
particles with opposite velocities at the center. The phase 
space structure at this time is presented in Fig O we no- 
tice that near the center we have two caustics. These two 
caustics are very close to the center and according to the 



theory of singularities for a one dimensional density, we ex- 
pect that their projected density will go like x~z , where x 
is a coordinate centered on the caustic. This is just what we 
obtain in Fig [4] The lower fold (red points in Fig [3J has a 
projected density close to x~? near the center. This density 
starting from the caustics extends at distances significantly 
larger than the separation between the caustic and the cen- 
ter, which means that in a large portion of this regime x ~ x 
and that as a consequence the density goes like x~?. The 
i~2 regime assume that the phase space density along the 
fold is constant, but due to the high non-linearity at the 
crossing the stretching of the phase space is large which im- 
plies that the density along the fold is nearly constant in a 
significant domain. This density associated with the caustic 
initiate the x~? asymptotic regime at the center and since 
it is dominant for /3o > — 5 it is prone to extend to larger 
values of x. At later times the system forms a large number 
of folds, with the persistence of the initial caustic induced 
x~? regime at the center. The x~ 5 central regime is also 
present for other initial conditions, for instance, a Gaussian, 
2 Gaussian,.., and many other, but the difference is that 
even if this regime is present very close to the center we find 
also a rather brutal transition to a tiny area with nearly 
constant density. This transition occurs at the very center 
of the system, and the size of this constant region depends 
on the initial conditions. This constant region at the center 
is probably a remain of the nearly flat density present in 
some initial conditions. The transition from the x~2 regime 
to the constant density regime signal also a break of the self 

similar regime. In any case the existence of this small con- 

_ i_ 

stant region affect the extent of the x 2 regime which is still 
present and most likely for the same reasons that for power 
law initial conditions. For a further investigation of this is- 
sue it is interesting to point the contribution of Colombi & 
Touma (2012) on this particular topic. 

9.2 Self similarity in the central region. 

It is important to note that in the central region of all nu- 
merical simulations self-similarity was observed (see Fig [8] 
for an illustration). The exponent — i of the density implies 
a.2 = 3, which means that in the x dimension the simulation 
should scale like j% and in the v dimension like . Using 
these scalings the folds at different times should be at the 
same location, which could be verified with good accuracy. 
This property definitely proves that the x~ 2 corresponds to 
a self similar regime in phase space. Since analytic formula 
for the self-similar solution with power law density were de- 
rived in the former sections, it is interesting to check the 
consistency with the fold system. The fold system is illus- 
trated in Figs [5]and [6] For power law initial conditions, it is 
expected that all quantities will be power laws, and thus the 
inter-fold distance should be given by Eq[2I] The potential 
4>(x) was computed, as well as R = \Jv 2 + cj>, and finally the 
inter-fold distance dR was estimated. The results are pre- 
sented in Fig [JJ The exponent of the inter- fold distance as a 
function of R according to Eq[2l]is 7 = ^- + 1 and provided 
that Q!2 is deduced from the initial condition exponent /3o, 
we have, a?2 = — -p- — 1 and 7 = (5^+2 • Fig [7] shows that this 
theoretical expectation of 7 is in good agreement with the 
numerical estimations. It is interesting to note that we ex- 
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Figure 2. A power law is fitted to the initial density and to the 
final late time density obtained after a large number of orbits. 
The resulting exponents are plotted as a function of each other 
for various initial conditions. For an initial exponent /3o < — \ the 
final exponent is identical, but for /3q > — \ a break occur and the 
final exponent saturates to —A, see text for more explanations. 
Note that in the range < /3o < — 1 the exponent was evaluated 
by fitting a power law to the force, and making the necessary 
conversion. 

pect that the theoretical expectation of 7 is going to match 
the data for /3o < — |, but for /3o > — \ we would expect 
/3 = — I and not /? = /3o. This is not what is observed, and 
we must conclude that despite the x _5 law at the center the 
structure of the fold system is preserved for /3q < — i. Other 
type of initial conditions imply different inter-fold distance, 
but this is always a small correction from the dR oc W law. 
To conclude it is interesting to note that all solutions for a 
large class of initial conditions have power law densities with 
exponents close to — i at the center, and are all self-similar 
with the same similarity index at the center. This definitely 
indicates that the degree of freedom (Q in Eq[l4]or H2 in 
Eq I15|) is important since it allows some flexibility to ac- 
commodate the initial conditions. All solutions have similar 
density at the center but have different phase space distri- 
butions while at the same time they remain self-similar. 



Figure 3. The phase space structure in the central region of the 
system just after the first crossings of the particles. This simula- 
tion contains a total of 4 X 10 5 particles. The initial conditions are 
dynamically cold with projected density profile x~3 and Gaus- 
sian velocity dispersion. 




10 PROPERTIES OF STATIONARY 

SELF-SIMILAR SOLUTIONS IN MULTI 
DIMENSIONS. 

The collapse of a dark matter halo in cosmological condi- 
tions is a very different problem from the ID collapse in- 
vestigated in the former section. The relation between the 
three dimensional collapse and self similarity was already 
explored by Lithwick & Dalai (2011). Here self similarity 
is imposed by the cosmological infall, the self similarity is 
not a consequence of a singularity at the center or of power 



Figure 4. The projected density associated with the first caus- 
tic. The first caustic corresponds to the red points in Fig [3] To 
demonstrate the relation between the power law projected den- 
sity and the projected density associated with the caustic, the x 
coordinate system was re-centered on the caustic position by in- 
troducing a shift in position xq. The red- line represents a power 
law with exponent — ~. 
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-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 

X 

Figure 5. The late time evolution of the phase space structure for 
the cold system. The initial conditions correspond to the density 
profile x~ 3 already presented in Fig [3] In this figure the most 
visible folds are the first outer folds, see Fig [6] for a more detailed 
view of the inner fold system. 




-4 -5 -2 -1 
ln(X) 

Figure 6. This plot presents a detailed view of the upper quad- 
rant (x > 0, v > 0) of the phase space diagram presented in Fig 
\E\ Logarithmic coordinates were introduced for both axis to allow 
a better view of the inner space folds. 

law initial conditions. In cosmology the initial stage of the 
halo formation can be represented with the accretion of mat- 
ter on a seed mass. This process implies that the turn-over 
radius of the in-falling shells of matter is a power law of 
time (Gunn 1977, Bertschinger, 1985). The relevant simi- 
larity exponent in the notations of this paper is «2 = — §■ 




nitial exponent 



Figure 7. The power law exponent 7 of the inter-fold distance 
dR. The exponent 7 was reconstructed by fitting the power law 
dR = R 1 to the fold system in phase space. The variable R is 
defined as: R = \/v 2 + <f>. The continuous line corresponds to the 
theoretical expectation (see Eq 121 H . 



I .5 




-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 

X 

Figure 8. The fold system at the center at two different times 
(black and blue). The initial conditions corresponds to the density 
profile x~ 3 already presented in Fig \3\ The blue folds are re- 
scaled using ct2 = 3. The time-scale factor is equal to the ratio 
of the time from the beginning of the simulation for each fold 
system. 
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The self similarity of the solution is enforced by the cos- 
mological infall, which implies that the final solution should 
reflect the similarity of the infall (Taylor & Navarro 2001). 
Thus we have good reasons to assume that the final out- 
come of the collapse of a dark matter halo will correspond 
to a self-similar solution of the Vlasov-Poisson system. Since 
dark matter is cold the solution will correspond to the mul- 
tiple folds of the initial surface in phase space. However, an 
equilibrium or a near equilibrium situation is a very differ- 
ent problem for a smooth continuous system in phase-space 
and for the folding of a cold system in phase space. Cold 
systems forms spiral in phase space and even in their late 
evolution they will be still evolving and will not converge 
to a smooth continuous equilibrium system. The smooth- 
ing of the spiral will produce a smooth phase-space density, 
and that smoothed density will be nearly stationary at late 
times. However in general the problem is that the smoothed 
density is not necessarily auto-similar. The smoothed den- 
sity will correspond to the smooth auto-similar solution only 
in some special cases. This problem is due to the fact that 
the smoothed phase space density does not satisfy the or- 
dinary Vlasov equation but a different differential equation. 
Since we cannot relate the cold self-similar solutions to the 
smooth self-similar solutions we cannot benefit from the fact 
that the smooth equilibrium solutions are simple analytic 
functions (Lancellotti & Kiessling, 2000). However it is clear 
that at late time the cold self similar solution should have a 
nearly stationary smoothed density in phase space and that 
it will correspond to a given sub-set of self similar solutions. 
The final problem is thus to identify the properties of this 
sub-set of solutions. A key point in the solution of this issue 
is to observe that even if the cold solution is discontinu- 
ous in phase space, the probability distribution of the phase 
space density P(f) could be similar to that of a smooth so- 
lution Ps(f). Actually if we consider the local surfaces of 
folds in phase space within some small solid angle, the inte- 
grated density between two local surfaces would correspond 
to the local surface multiplied by the thickness of the fold. 
The integration and averaging of the density between two 
folds corresponds to the smoothing operation. If the fold 
thickness is proportional to the distance between the folds, 
then the smoothed density will be proportional to the cold 
un-smoothed density, and its occurrence will also be pro- 
portional to the cold occurrence since it is proportional to 
the ratio between the thickness and the inter-fold distance. 
This reasoning shows that the smoothed probability density 
of / can be deduced from the un-smoothed density by the 
introduction of constant scale factors (see Fig [9] for an illus- 
tration) . As a consequence the self-similar properties of the 
cold probability distribution will be also valid for the smooth 
probability distribution. Note that / is made of many dif- 
ferent contours, but since the scaling properties between P 
and Ps are valid for each contour, the transfer of self simi- 
larity between P and Ps remain valid in presence of many 
different contours. The property of proportionality between 
the thickness and the inter fold distance was demonstrated 
in Sec. 6, but is also required for the existence of the solu- 
tion. If the thickness exceeds the inter fold distance then a 
crossing would occur near the origin, while if the opposite 
behavior happened the crossing would happen on larger dis- 
tances. The only solution is thus that the fold thickness is 
proportional to the fold inter distance. As a consequence it 



is clear that the smooth probability distribution inherit the 
self similar properties of the cold solution, and that power- 
law integrated quantities estimated from Ps(f) should have 
an auto-similar behavior. Note that a priori the equivalence 
between P(f) and Ps(f) is valid on any sub domain of the 
phase space, provided that this domain is large enough to 
define the statistics of /. In particular all quantities like: 

Qn(x )= j P S (/(x,v))/(x,v)"df 

with : x = xo and : < v < oo 
have to be self similar since 

Q°(x ) = j P(/(x,v))/(x,v)"d/ 

with : x = xo and : < v < oo 

Is self-similar and Ps inherits the self similarity of P. As 
a consequence the smooth self-similar Q n function should 
write as a power law of time multiplied with a function of 
X2. A direct estimate of the mean value of / Qi is: 

In the same spirit one may approximate Q n : 



The ratio of Q n quantities is also self-similar: 




Since Q„ and R„ m are smooth continuous functions and self 
similar: 

Q n (x)=t na °hi(x 2 ) 

and 

P„ m (x) = t ( "- m)ao Mx 2 ) 

Stationary solutions for Q n and R nm implies that /ti and I12 
are power law's, thus: 

Qn(r) oc t na ° [^-] 71 = r"' 1 ; 71 = n^- = -n^ 

Similarly: 

^W«t ( "" mH [^f 2 =^ ; T> = -(n-m)y 

(Note that we assume isotropic systems, thus x was replaced 
with the distance modulus r.) 

The calculation of integrals involving higher order power 
of f (f n , n > 1) requires the knowledge of the spatial 
distribution of /. The reconstruction of / from numerical 
discrete numerical data is a difficult and uneasy task (See 
Arad, Dekel, & Klypin, 2004). A more convenient way of 
estimating the higher order terms is to use the accurate 
analytic approximations of / that were developed by 
Wojtak, etal (2008). All integrals were performed using Eq 
29 from Wojtak, etal (2008). The results are presented in 
Figs [TO] and 1111 For convenience an inverse power law's of 
higher order quantities was taken to reduce all self similar 
exponents to the — theoretical self-similar expectation. 
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All quantities are nearly linear in the log-log diagram, and 
all slopes in this diagram are very close to — The only 
possible explanation is that if all moments of the smooth 
probability distribution Ps(f) are self similar is that Ps(f) 
itself is self similar. As a consequence the universal behavior 
observed for CDM haloes in numerical simulations reflects 
the self similarity of Ps(f). 



11 CONCLUSION AND PERSPECTIVES. 

The convergence of cold initial conditions with initial 
density profile not steeper that x~ 5 to a universal power 
law density with exponent — | (the Binney conjecture, BC) 
is an interesting and powerful example of universality. The 
one dimensional solutions are quite simple, and the BC is a 
very useful and efficient tool to explore self similarity. The 
conjecture implies that a large class of self similar solutions 
forms at the center of collapsed systems. The analytic solu- 
tions developed in this paper help to understand the general 
properties of the solutions, and the basic mechanism at the 
origin of the conjecture, which is related to a singularity has 
been identified. However the full richness and complexity 
of the ID self similar solutions have yet to be explored, 
and the BC provides a very interesting method to generate 
a large class of self similar solutions. The final self similar 
solution will depend on the initial conditions, and it is 
clear that more complicated and exotic initial conditions 
have to be explored. The problem of the 3D cold collapse 
in a cosmological context, is somewhat different since the 
similarity is not due to a singularity but is enforced by the 
nature of the infall. In 3D simulations of dark matter haloes 
the pseudo phase space density is an universal power 
law. It was proved that this similarity is not the only one, 
and that other moments of the probability distribution of 
/ have the same universal properties. The convergence to a 
power law of all these quantities is due to the self similarity 
of the probability distribution of the smoothed /, which is 
inherited from the self similarity of the original distribution. 
Interestingly it was shown by Dehnen & McLaughlin (2005) 
that the universal power law implied universal density 
profiles, not exactly NFW profiles, but quite close. It is 
interesting to point that the NFW profile itself may be 
understood usimg simple physics (Dalai, N., Lithwick, Y., 
Kuhlen, M., 2010). As a consequence one may conclude 
that all the universality properties of the dark matter haloes 
are related the self similarity of the probability distribution 
of the smoothed /, Ps(f). A final point is that as expected, 
the self similar properties of Ps(f) are also present for 
the one dimensional solutions presented in this paper, and 
that higher order moment of Ps(f) are also power laws, 
which is a fundamental property of self similar solutions. 
A fundamental and natural perspective is to extend this 
work to real galaxies where not only dark matter is present 
but also baryonic matter. It is remarkable that the rotation 
curves of galaxies show also some universal similarity 
properties (Donato etal 2009, Salucci etal 2007.) 

The author would like to thank Scott Tremaine for 
support during his stay at the Institute for advanced 



Figure 9. The cold un-smoothed density is represented in blue 
and the local value of the smoothed density is represented in red 
in this simple ID analogy. The cold density is at the position of 
the fold and its thickness is proportional to the inter-fold dis- 
tance. The smooth density corresponds to the averaging of the 
cold density between two consecutive folds. Since the thickness 
and the inter-fold distance are proportional, the construction of 
the smooth density is equivalent to a simple re-scaling of the cold 
density. Obviously the cold and smoothed density cannot be pro- 
portional in all points of phase space, since the cold density is zero 
in some place, and the smooth density is neither equal to zero be- 
tween folds. However the distribution of the values of the smooth 
and cold density differs only by scale factors, which implies that 
the self similar (scaling) properties of the cold distribution of / 
values are preserved in the smooth density distribution, or prob- 
ability distribution of the smooth /. 
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Figure 10. The following quantities black, 



ff 2 d 3 v 



0.0 



blue, 



If 3 d 3 v 



red, are represented by points 

(cross). The points were calculated using the Wojtak etal (2008) 
model. A straight line with slope — ^ is over-plotted over each 



studies in 2011. The author thanks Stephane Colombi, 
Jean Philippe Beaulieu, and Jihad Touma for interesting 
comments. 
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